function [xi0,xi1,xis,xim,xip0,xip1,xip,phi0,phi1,phis,phip0,phip1,phip,Cs,Cp,u,Fs]=Solve_model_E(Paras)
%%%%%��������%%%%%%%%%%
 global xb0 xb1  phi0 phi1 phis xi0 xi1 xis xip xim phip0 phip1 phip
     r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    P=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
    c=P*r/(1-taui);
    
    psi=mu1-mu0;
    mu=mu0+psi*pi0;
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
       
    beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
    
    
    
    %%%%�Ʋ�ˮƽ%%%%%
    xb0=(beta2/(beta2-1))*((r-mu0)/r)*c;
    xb1=(beta4/(beta4-1))*((r-mu1)/r)*c;
    
    
    %%%%%Ͷ��ˮƽ%%%%%%%%%%%%
     y0=0.8;
     y0=fsolve(@(P)FB(P,Paras),y0);
 
     xi0=roundn(y0,-8); 
     
     
     %%%%%%%������%%%%%%%%%%%%%%
     E0=(1-taue)*(xi0/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(xi0/xb0)^beta2);
     D0=(1-taui)*c/r+((1-alpha)*(1-taue)*xb0/(r-mu0)-(1-taui)*c/r)*(xi0/xb0)^beta2;
     
     phi0=((1-taue)*c/r-(1-taue)*D0/(1-taui))/E0;
     
     
   %%%%%Ͷ��ˮƽ%%%%%%%%%%%%
     y0=0.8;
     y0=fsolve(@(P)SB(P,Paras),y0);
 
     xi1=roundn(y0,-8); 
     
     
     %%%%%%%������%%%%%%%%%%%%%%
     E1=(1-taue)*(xi1/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xi1/xb1)^beta4);
     D1=(1-taui)*c/r+((1-alpha)*(1-taue)*xb1/(r-mu1)-(1-taui)*c/r)*(xi1/xb1)^beta4;
     
     phi1=((1-taue)*c/r-(1-taue)*D1/(1-taui))/E1;
     
     Fg=((1-phi1)*E1-(I-(1-taui)*c/r))*(x/xi1)^beta3;
     
      %%%%%����Ͷ��ˮƽ%%%%%%%%%%%%
      y0=1;
      y0=fsolve(@(P)SEB0(P,Paras),y0);
 
     xis=roundn(y0,-8);  
     
    
     
          %%%%%%%������%%%%%%%%%%%%%%
     e1=(1-taue)*(xis/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xis/xb1)^beta4);
     d1=(1-taui)*c/r+((1-alpha)*(1-taue)*xb1/(r-mu1)-(1-taui)*c/r)*(xis/xb1)^beta4;
     
     phis=((1-taue)*c/r-(1-taue)*d1/(1-taui))/e1;
     
     Fs=((1-phis)*e1-(I-(1-taui)*c/r))*(x/xis)^beta3;
     

     
       
     
     
               %%%%%badtypepoolingͶ��ˮƽ%%%%%%%%%%%%
      y0=1;
      y0=fsolve(@(P)PB0(P,Paras),y0);
 
     xip0=roundn(y0,-8);  
    
     
                    %%%%%goodtypepoolingͶ��ˮƽ%%%%%%%%%%%%
      y0=1;
      y0=fsolve(@(P)PB1(P,Paras),y0);
 
     xip1=roundn(y0,-8);  
     
                    %%%%%%%%%%%%%%�������poolingͶ��ˮƽ 
     xig=xis:0.0001:xi0;
     n=length(xig);
     V=zeros(n,1);
     phip=zeros(n,1);
     pi=zeros(n,1);
     u=zeros(n,1);
     
     for j=1:n
     
         
       xip=xig(j);
      
      pi(j)=1/(1+(1-pi0)*(x/xip)^(psi/sigma^2)/pi0);
     % pi0(j)=pi(j);
         
     %%%%%%%%%%%�Ʋ�ˮƽ%%%%%%%%%%%%%%
  %  xb(j)=(((1-pi(j))*beta2+pi(j)*beta4)/((1-pi(j))*(beta2-1)/(r-mu0)+pi(j)*(beta4-1)/(r-mu1)))*(c/r);
    
         %%%%%%%������%%%%%%%%%%%%%%
     E0(j)=(1-taue)*(xip/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(xip/xb0)^beta2);
     D0(j)=(1-taui)*c/r+((1-alpha)*(1-taue)*xb0/(r-mu0)-(1-taui)*c/r)*(xip/xb0)^beta2;
     
     E1(j)=(1-taue)*(xip/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip/xb1)^beta4);
     D1(j)=(1-taui)*c/r+((1-alpha)*(1-taue)*xb1/(r-mu1)-(1-taui)*c/r)*(xip/xb1)^beta4;
     
     phip(j)=((1-taue)*c/r-(1-taue)*((1-pi(j))*D0(j)+pi(j)*D1(j))/(1-taui))/((1-pi(j))*E0(j)+pi(j)*E1(j));
    
    
    
    %%%%%%%��Ȩ��ֵ%%%%%%
    Em(j)=(1-taue)*(xip/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip/xb1)^beta4);
   % x=0.5; 
    
    V(j)=((1-phip(j))*Em(j)-(I-(1-taui)*c/r))*(x/xip)^beta3;   
         
    u(j)=V(j);     
     end
     [I]=find(V==max(V));
  
     xip=xig(I);
     phip=phip(I);
     Fp=V(I);
     pi=pi(I);
     u=u(I);
     
     
        %%%%%��СͶ��ˮƽ%%%%%%%%%%%%
      y0=1;
      y0=fsolve(@(P)SEB1(P,Paras),y0);
 
     xim=roundn(y0,-8);  
     
    %%%%%%%%%�ɱ�%%%%%%%%%%
    Cs=(Fg-Fs)/Fg;
    Cp=(Fg-Fp)/Fg;
   
     
     
     
     
     
     


end

function F=FB(P,Paras)
 xi0=P;
%%%%%��������%%%%%%%%%%
 global xb0 
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    P=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
    c=P*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    
    

     F=((taue-taui)*c/r-I)*beta1+(beta1-1)*(1-taue)*xi0/(r-mu0)+((beta1-beta2)*((1-taue)/(1-taui))*((1-taue)*(1-alpha)-(1-taui))*xb0/(r-mu0))*(xi0/xb0)^beta2;


end

function F=SB(P,Paras)
 xi1=P;
%%%%%��������%%%%%%%%%%
 global xb1 
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    P=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
    c=P*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
    

     F=((taue-taui)*c/r-I)*beta3+(beta3-1)*(1-taue)*xi1/(r-mu1)+((beta3-beta4)*((1-taue)/(1-taui))*((1-taue)*(1-alpha)-(1-taui))*xb1/(r-mu1))*(xi1/xb1)^beta4;


end

function F=SEB0(P,Paras)
 xis=P;
%%%%%��������%%%%%%%%%%
 global xb0 xb1 phi0 phi1 xi0 xi1
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    P=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
    c=P*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
      beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    %%%%%%%��Ȩ��ֵ%%%%%%
   e1=(1-taue)*(xis/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xis/xb1)^beta4);
    d1=(1-taui)*c/r+((1-alpha)*(1-taue)*xb1/(r-mu1)-(1-taui)*c/r)*(xis/xb1)^beta4;
     
     phis=((1-taue)*c/r-(1-taue)*d1/(1-taui))/e1;
     
     Es=(1-taue)*(xis/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(xis/xb0)^beta2);
    E0=(1-taue)*(xi0/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(xi0/xb0)^beta2);
    

   F=(1-phis)*Es-(I-(1-taui)*c/r)-((1-phi0)*E0-(I-(1-taui)*c/r))*(xis/xi0)^beta1;


end

function F=SEB1(P,Paras)
 xim=P;
%%%%%��������%%%%%%%%%%
 global xb0 xb1 phi0 phi1 phip xi0 xi1 xip
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    P=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
    c=P*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
    %%%%%%%��Ȩ��ֵ%%%%%%
    Em=(1-taue)*(xim/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xim/xb1)^beta4);
    
     E1=(1-taue)*(xip/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip/xb1)^beta4);
 
    
    
    

   F=(1-phi1)*Em-(I-(1-taui)*c/r)-((1-phip)*E1-(I-(1-taui)*c/r))*(xim/xi0)^beta3;


 


end

function F=PB0(P,Paras)
 xip0=P;
%%%%%��������%%%%%%%%%%
 global xb0 xb1 phi0 phi1 xi0 xi1 phip0
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    P=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
    c=P*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
       
    beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
%     %%%%%%%%%%%�Ʋ�ˮƽ%%%%%%%%%%%%%%
%     xb=(((1-pi0)*beta2+pi0*beta4)/((1-pi0)*(beta2-1)/(r-mu0)+pi0*(beta4-1)/(r-mu1)))*(c/r);
%     
         %%%%%%%������%%%%%%%%%%%%%%
     E0=(1-taue)*(xip0/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(xip0/xb0)^beta2);
     D0=(1-taui)*c/r+((1-alpha)*(1-taue)*xb0/(r-mu0)-(1-taui)*c/r)*(xip0/xb0)^beta2;
     
     E1=(1-taue)*(xip0/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip0/xb1)^beta4);
     D1=(1-taui)*c/r+((1-alpha)*(1-taue)*xb1/(r-mu1)-(1-taui)*c/r)*(xip0/xb1)^beta4;
     
     phip0=((1-taue)*c/r-(1-taue)*((1-pi0)*D0+pi0*D1)/(1-taui))/((1-pi0)*E0+pi0*E1);
    
    
    
    %%%%%%%��Ȩ��ֵ%%%%%%
    Es=(1-taue)*(xip0/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(xip0/xb0)^beta2);
    E0=(1-taue)*(xi0/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(xi0/xb0)^beta2);
    

   F=(1-phip0)*Es-(I-(1-taui)*c/r)-((1-phi0)*E0-(I-(1-taui)*c/r))*(xip0/xi0)^beta1;


end

function F=PB1(P,Paras)
 xip1=P;
%%%%%��������%%%%%%%%%%
 global xb0 xb1 phi0 phi1 phis xi0 xi1 xis xim phip1
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    P=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
    c=P*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
       
    beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
%         %%%%%%%%%%%�Ʋ�ˮƽ%%%%%%%%%%%%%%
%     xb=(((1-pi0)*beta2+pi0*beta4)/((1-pi0)*(beta2-1)/(r-mu0)+pi0*(beta4-1)/(r-mu1)))*(c/r);
    
         %%%%%%%������%%%%%%%%%%%%%%
     E0=(1-taue)*(xip1/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(xip1/xb0)^beta2);
     D0=(1-taui)*c/r+((1-alpha)*(1-taue)*xb0/(r-mu0)-(1-taui)*c/r)*(xip1/xb0)^beta2;
     
     E1=(1-taue)*(xip1/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip1/xb1)^beta4);
     D1=(1-taui)*c/r+((1-alpha)*(1-taue)*xb1/(r-mu1)-(1-taui)*c/r)*(xip1/xb1)^beta4;
     
     phip1=((1-taue)*c/r-(1-taue)*((1-pi0)*D0+pi0*D1)/(1-taui))/((1-pi0)*E0+pi0*E1);
    
    
    
    %%%%%%%��Ȩ��ֵ%%%%%%
    Em=(1-taue)*(xip1/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip1/xb1)^beta4);
    E1=(1-taue)*(xis/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xis/xb1)^beta4);
    

   F=((1-phip1)*Em-(I-(1-taui)*c/r))*(1/xip1)^beta3-((1-phis)*E1-(I-(1-taui)*c/r))*(1/xis)^beta3;

end



